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ABSTRACT 



> 

' We describe a new hybrid iV-body/hydrodynamical code based on the particle- mesh 

m ' 

CO 

■ related to the evolution of large-scale structure, galaxy clusters, and individual galaxies. 

The code, named COSMOS, possesses several new features which distinguish it from 
, other PM-PPM codes. In particular, to solve the Poisson equation we have written a 

new multigrid solver which can determine the gravitational potential of isolated matter 
distributions and which properly takes into account the finite-volume discretization re- 
quired by PPM. All components of the code are constructed to work with a nonuniform 

S-H ' 

mesh, preserving second-order spatial differences. The PPM code uses vacuum bound- 
ary conditions for isolated problems, preventing inflows when appropriate. The PM 
code uses a second-order variable-timestep time integration scheme. Radiative cooling 
and cosmological expansion terms are included. COSMOS has been implemented for 
parallel computers using the Parallel Virtual Machine (PVM) library, and it features 
a modular design which simplifies the addition of new physics and the configuration of 
the code for different types of problems. We discuss the equations solved by COSMOS 
and describe the algorithms used, with emphasis on these features. We also discuss the 
results of tests we have performed to establish that COSMOS works and to determine 
its range of validity. 

Subject headings: Methods: numerical — hydrodynamics — dark matter — gravitation 
— large-scale structure of universe — intergalactic medium 
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1. Introduction 

The study of many problems of interest in astrophysics today benefits from or even requires 
the use of three-dimensional hydrodynamical or iV-body simulations because the complex physical 
processes or the lack of symmetry involved prohibit analytic or reduced-dimensionality approaches. 
These problems span the entire range of astrophysical length scales, from black hole collisions and 
core-collapse supernovae to clusters of galaxies and the Lyman a forest. The advent over the past 
two decades of modern shock-capturing hydrodynamical algorithms and fast iV-body methods, 
together with rapid advances in available computing power, has made such simulations feasible. 

In these simulations one must make a distinction between collisional matter, in which the 
mean free path between particle collisions is much smaller than all length scales of interest, and 
collisionless matter, in which particles may free-stream between collisions on length scales important 
to the simulation. For collisional matter, a fluid description incorporating a pressure field and a 
single-valued bulk velocity field is appropriate; shock waves can convert energy irreversibly from 
bulk kinetic form to internal form. For collisionless matter, the bulk velocity can take on multiple 
values at each point in space, and shock waves are not possible. 

For simulations of galaxies, clusters of galaxies, and large-scale structure it is important to 
track both kinds of matter. On the length and time scales of interest in these simulations, baryonic 
matter (the intracluster or intergalactic medium) behaves as a collisional fluid. Although locally 
this gas may be partially or fully ionized, with collisions between ions and electrons occurring by 
means of randomly oriented magnetic fields, on resolvable length scales we may assume charge 
neutrality and ignore local deviations from fluid behavior (e.g., plasma instabilities). However, 
because we do not yet know the identity of the dark matter, we adopt the simplest hypothesis and 
assume that it is collisionless. The gravitational stability of apparently bound systems (galaxies 
and clusters of galaxies) in which baryons provide insufficient mass is still our best evidence for the 
existence of dark matter. Thus the simplest hypothesis also treats the interaction between these 
types of matter as purely gravitational. Cosmological simulations must track the evolution of both 
matter constituents in their mutual gravitational field. 

Several codes now exist to simulate the combined evolution of collisional and collisionless matter 
in a cosmological setting. All treat the dark matter using superparticles which sample the particle 
distribution. These iV-body codes differ in the methods they use to compute interparticle forces 
and to solve the hydrodynamical equations for the gas. Methods for computing interparticle forces 
include direct summation (PP, or particle-particle), tree algorithms (Hernquist 1987), particle-mesh 
methods (PM: Hockney & Eastwood 1988), and combinations of direct summation and particle- 
mesh (P 3 M: Efstathiou & Eastwood 1981). The tree and particle-mesh algorithms trade off force 
accuracy or spatial resolution for speed in comparison with direct summation, while variants of 
P 3 M attempt to overcome the resolution limitations of PM techniques while retaining their speed. 

The hydrodynamic solvers currently in use with particle codes include variants of smoothed- 
particle hydrodynamics (SPH: Gingold & Monaghan 1977; Lucy 1977) and grid-based methods 
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such as the piecewise-parabolic method (PPM: Colella & Woodward 1984) and the total-variation- 
diminishing method (TVD: Harten 1983). SPH algorithms descended from astrophysical iV-body 
techniques and use particles to approximate the behavior of the gas, treating gas particles as 
moving interpolation centers for quantities such as the gas pressure. Typically SPH codes achieve 
good spatial resolution in high-density regions but handle shocks and low-density regions poorly. 
Examples of cosmological hydrodynamics codes based on SPH include those of Evrard (1988), 
Hernquist & Katz (1989), Navarro & White (1993), Couchman, Thomas, & Pierce (1995), Steinmetz 
(1996), and Owen et al. (1998). Grid-based methods, which are used widely in other areas of 
physics and in engineering, suffer from more limited resolution, but they handle high-density and 
low-density regions equally well, and with modern algorithms they handle shocks extremely well. 
Moreover, it is more straightforward to add extra gas physics to grid-based codes. Examples of 
grid-based cosmological hydrodynamics codes (excluding those based on PPM, which are discussed 
below) include the centered-difference code of Cen (1992), the TVD code of Ryu et al. (1993), 
the softened Lagrangian code of Gnedin (1995), and the moving-mesh TVD code of Pen (1998). 
Comparisons of various cosmological particle- and grid-based codes have been performed by Kang 
et al. (1994) and Frenk et al. (1999). 

The combination of PM for dark matter and PPM for gas has proven to be an especially useful 
method for cosmological simulations. Accurate shock handling and straightforward implementation 
of subgrid physics argue for the use of a grid-based scheme for the gas. When using a grid for the 
gas, the most efficient means for obtaining the forces on the dark matter particles is the particle- 
mesh method, because the gas and dark matter both make use of the potential defined on the grid. 
Also, because the gas grid limits the spatial resolution, the greater dynamic range of a tree code 
or direct summation solver for the dark matter is not as useful. In recent years several groups 
have independently developed PM-PPM codes to study large-scale structure and galaxy cluster 
formation (Bryan et al. 1995; Sornborger et al. 1996; Gheller, Pantano, & Moscardini 1998; Quilis, 
Ibahez, & Saez 1996, 98) and cluster evolution (Roettiger, Stone, & Mushotzky 1997). These 
codes differ in a number of ways. The Bryan et al., Sornborger et al., and Roettiger et al. codes 
use the Lagrange-plus-remap formulation of PPM, whereas the Gheller et al. code uses the direct 
Eulerian formulation. Little difference in numerical accuracy between the two formulations has 
been observed, but the Lagrange-plus-remap formulation generalizes somewhat more readily to 
moving or adaptive grids. The Bryan et al. code adds certain ad hoc modifications to the basic 
PPM scheme to improve the resolution of narrow density peaks, and more recently this code has 
been extended to include adaptive mesh refinement (AMR: Norman &: Bryan 1999). The Quilis 
et al. code makes use of a linearized Riemann solver due to Roe (1981), whereas the other codes 
use variants of the iterative Riemann solver due to van Leer (1979). All of the codes except that 
of Roettiger et al. use a Green's function technique (via the Fast Fourier Transform, or FFT) to 
solve for the combined gravitational potential of the gas and dark matter. The Roettiger et al. 
code uses a conjugate gradient solver (Hockney & Eastwood 1988) for the Poisson equation. To 
our knowledge, all of the codes use cloud-in-cell (CIC) weighting together with a variable-timestep 
time integration scheme for the particle-mesh code. 
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We have developed a new PM-PPM code, called COSMOS, for simulations of large-scale 
structure formation and galaxy cluster evolution. Like Gheller et al. (1998), we use the direct 
Eulerian formulation of PPM, and like all of the above codes we use CIC weighting and a variable 
timestep in our PM code. However, our code uses a full multigrid solver rather than an FFT to 
obtain the gravitational potential, enabling us to handle nonuniform grids with isolated boundary 
conditions. Also, we have implemented a nonideal equation of state following Colella & Glaz 
(1985) (of use, for example, in Ly a cloud and supernova simulations), and we implement cooling 
and multifluid source terms, as well as cosmological expansion terms, with a method similar to that 
used by Kritsuk, Bohringer, & Miiller (1998). The hydrodynamical and gravitational components 
of COSMOS have been used previously by Ricker (1998) to study purely hydrodynamic cluster 
evolution, and the full code has been used by Ricker & Sarazin (2000) to simulate cluster mergers 
including gas and dark matter. 

In this paper we describe the COSMOS code and the results of our validation tests. Section 
2 discusses the physics included in the code and the corresponding equations. In Section 3 we 
discuss the discretization scheme used for the grid-based components of the code, the algorithms 
we use, and the modifications we have made to the standard algorithms to suit the requirements 
of cosmological problems. Section 4 describes the test problems we have used to validate the code. 
Finally, Section 5 summarizes the paper. 



2. Equations 

In this section we describe the definitions and equations used in the COSMOS code. All 
calculations take place in a three-dimensional computational volume in which positions are specified 
using comoving coordinates x = r/a, where r is a proper position vector and a(t) is the time- 
dependent cosmological scale factor. Both dark matter and gas peculiar velocities v are considered 
to be comoving; thus v = x for a particle at location x. We define the comoving gas density and 
pressure as 

Pg = a3 Pg (!) 
p = ap , 

where p g and p are the proper gas density and pressure, respectively. (We use the subscripts g and 
dm to distinguish quantities which apply separately to both the gas and dark matter.) Note that 
our definition of p requires that the comoving internal energy density be given by 

p g e = ap g e , (2) 

where e is the proper specific internal energy of the gas. For an ideal-gas equation of state with 
ratio of specific heats 7 and average atomic mass fi, we have 

Pg e = 7 = / 7\ , ( 3 

7-1 (7-l)A» 
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where 

f 

T=- (4) 

cr 

is the comoving temperature, and T is the proper temperature. We also define the comoving 
potential 4> as the solution at each simulation timestep to Poisson's equation, 

V 2 6 = 



a? 



(Pg + Pdm) - (Pg + Pdm) • (5) 



Here pdm is the comoving density in dark matter particles, and the bar indicates a spatial average. 
The comoving total energy density of the gas is defined as 

PgE = p g e + lp g vl + p g <f> , (6) 

where v s = |v g | is the comoving speed of the gas. Finally, for certain problems we make use of an 
equilibrium cooling function, which we denote by A(p g ,T). For example, cooling via optically thin 
thermal bremsstrahlung emission, assuming a fully ionized gas, would require A oc p g 2 T 1 / 2 . 

We use these definitions to transform the inviscid Eulerian equations of hydrodynamics into 
a simple comoving form which, in a static universe (a = 1), reduces again to the standard equa- 
tions with comoving quantities equivalent to proper ones. The comoving gas equations solved by 
COSMOS are as follows: 

d 

-p g + V • [p g v g ] = (7) 

— (p g v g ) + V • [pgv g v g ] + Vp + 2% gV g + p g V4> = (8) 

- (p g E) + V • [( Pg E + p) v g ] + - a [(3 7 - l)p g e + 2 Pg v 2 g ] - Pg ^ + A = (9) 
d ex 

— (p g e)+V - [(pge + p)vg]-Vg-Vp+^(3 7 -l)p g e + A = (10) 

The equations are written in explicit conservation form. With the finite-volume discretization 
used in COSMOS (Section 3.1.1), we can difference these equations (except for the internal energy 
equation) in such a way that errors in the advection terms do not affect global conservation of 
mass, momentum, and energy. 

We determine which energy equation to use in updating the pressure on the basis of the 
magnitude of the internal energy density. During a simulation timestep we update the internal 
energy and pressure according to the internal energy equation (10) and the equation of state (3). 
We update the total energy using the total energy equation (9). At the end of the timestep, if in 
any zone the condition 

Pg e > 10" 2 max jipgv 2 , |p g 0|j (11) 

is satisfied, we reset the internal energy density using equation (6) with the updated values of p g E, 
p g , v g , and 4>, then use the equation of state to reset the pressure. If this criterion is not satisfied, 
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we instead reset the total energy density in that zone using equation (6). This procedure permits us 
to maintain total energy conservation in regions where the pressure makes a significant contribution 
to the total energy while avoiding large round-off errors in the pressure where it is small. Condition 
(11) is generally satisfied except for high Mach-number flows and in 'dust' calculations where the 
gas pressure offers little support against gravity. Similar methods based on energy (Bryan et al. 
1995) and entropy (Ryu et al. 1993) are used by other codes to handle such flows. 

Equations (7) through (10) suffice to describe the behavior of matter that is collisional on the 
length scales of interest to us. Examples include the gas in stellar atmospheres and in clusters of 
galaxies. However, for collisionless matter components, such as galaxies and dark matter in clusters, 
we must solve the equations of motion for individual particles. These particles are subject only to 
gravitational forces; hence a particle's comoving position Xd m and velocity Vd m evolve according to 

^ x dm 



dt 



Vdm (12) 



— — + 2-v dm = -V0, (13) 
dt a 

the second term on the left in the velocity equation representing the Hubble flow due to uniform 
expansion. 



3. Numerical methods 
3.1. Hydrodynamics 

3.1.1. Finite-volume discretization 

We solve the hydrodynamical equations described in Section 2 on a nonuniform, finite, Carte- 
sian grid containing N x x N y x N z cells. The center of the (i, j, k)-th cell (i = 1 . . . N x , j = 1 . . . Ny, 
k = l...Nz) is located at x^ =(xi,yj, Zk). The edges of this cell have x-coordinates a?i±i/2 = 
Xi ± T^Axi and y- and ^-coordinates that are similarly defined. 

We use a standard finite- volume discretization. A convenient way of expressing this discretiza- 
tion is as follows. We define a convolution operator which averages a variable over a cell volume: 

i f-x+Axi/2 ry+Ayj/2 f -z+Az k /2 

<*Vk\fW] = -T T T / / / dV/(x'). (14) 

£\Xi lAyj l\Z k Jx-Axi/2 Jy-Ayj/2 J z-Az k /2 

Here / is any scalar fluid variable (e.g. density or pressure). Now let 

fijk = ctij k [f(xijk)} ; (15) 
by substituting a Taylor expansion for /(x') about x^ we find that 

fijk = /(Xi jfc ) + 0(A 2 ) . (16) 
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Thus the fluid quantities manipulated by this finite-volume technique differ from those used in 
finite-difference methods by terms of order the square of the grid spacing. These cell averages are 
to be distinguished from cell-wall averages, which we define as follows for cell walls perpendicular 
to the x-axis: 

AyjAzkJy._Ay.f2 Jz k -Az k /2 
Similar definitions apply for the other coordinate directions. 



I ryj+^yj/ 2 r^k+^ k /2 
Si±\/2_k = -z T— / dy'dz' f(x i±1/2 ,y',z') . (17) 

A yj Az k J v .- Av ./ 2 Jzu-Az h /2 



If we write the fluid equations (7) - (10) with x',y',z' as the spatial variables, apply the 
convolution operator (14) to both sides, apply the divergence theorem, and then let x = Xj, y = yj, 
and z = Zk, we obtain a set of spatially discretized equations. As usual, the continuity equation 
becomes 

d 1 

di Ps ' ijk+ ~Ax~ [(Ps V g,'^)i+l/2,jk ~ (PgVg,x)i-l/2,jk\ + 

[(Ps v s,v)i,j+i/2,k - (.Pg v g,y)i,j-l/2,k\ + (18) 

^J" [{PgVg,z)i,j,k+l/2 - {PgVg,z)i,j,k-l/2\ = • 

The other fluid equations transform in a similar fashion. In this form the advection terms in 
equations (7) - (9) are conservatively differenced, because the same amount of each conserved 
quantity is subtracted from each cell as is added to its neighbors. The nonconservative advection 
term in the internal energy equation (10) does not significantly affect global energy conservation, 
because wherever the internal energy density makes a large contribution to the total it is reset 
using the total energy density, which is conservatively differenced. 

Because the hydrodynamical equations are hyperbolic, we solve them by integrating forward 
in time from a given initial state. We use an explicit forward time difference, so for stability we 
limit the timestep using the Courant-Friedrichs-Lewy (CFL) criterion (Roache 1998): 

Axi Ayj Az k 



Athydro = o min 

grid 



|^g,x,ijfc| ~b |Cs,jj'fc| l^g,j/,ijfc| "b |Cs,jj'fc| |^g,2,ijfc| |Cs,jj'fc| 

Here a is the 'CFL parameter', a number between and 1 (as close to 1 as possible for accuracy 
and for computational efficiency), and the minimization is taken over the entire computational grid. 
Additional restrictions on the timestep due to gravitational acceleration (At grav ), radiative cooling 
(Aicooi), and the iV-body code (Ai^m) are described in the following sections. We adopt a timestep 
At which is the minimum of all of these restrictions: 

At = min [At hydro , At grav , At coo i, At dm , . . .] . (20) 

To obtain the fully discretized equations we average the spatially discretized equations over 
the time interval (t n , t n+ \), where At n = t n+ \ — t n . The continuity equation then becomes 

At r , 



Pg^ijk ~ Pg,ijk ^ gUg ' x ^+ 1 / 2 .J'fe (P& V &> x )i-l/2,jk 
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At 



n 



At n 



(P& v g,y)i,j+i/2,k (Ps v s,y)i,j-i/2,k ( 21 ) 



(Pg v S,z)iJ,k+l/2 (Pg v g, z h,j,k-i/2 



where bars indicate averages over (t n ,t n+ i). This is an exact equation for Pg^^.; discretization 
error is introduced when we estimate the time-averaged fluxes ({PgV gjX ) i+1 / 2 jfe' e * c -) gi yen the cen_ 
averaged fluid variables at time t n , since the fluxes depend on cell-wall averages which must be 
determined through interpolation. In COSMOS the interpolation and flux computation steps are 
handled using the piecewise-parabolic method (PPM), which we discuss in the next section. 

To simplify the algorithm, we reduce this three-dimensional problem to a set of one-dimensional 
problems using standard operator splitting techniques. This allows us to use a one-dimensional 
hydrodynamics routine which operates on ID work arrays, into which and out of which we swap 
rows or columns from the full grid as necessary. In the simplest form of operator splitting, one 
advances every cell through At using just the x-derivatives, then takes the results of this operation 
and advances them again through At using just the ^-derivatives, finally repeating the procedure 
using just the z-derivatives. This method produces an effective three-dimensional operator which 
is accurate to O(At). We have found in spherical explosion tests that first-order splitting does not 
preserve rotational symmetry well. Therefore we use the symmetrized operator splitting method 
of Strang (1968), which yields 0(At 2 ) accuracy. For each successive pair of timesteps we perform 
first a sweep in x, then a sweep in y, then a sweep in z, each for a full timestep. We then repeat 
these operators in reverse order for the second timestep. 

Splitting techniques also generalize to the other operators in the problem. For the cosmological 
expansion terms we use the exact solution of 

^(Pgv g ) = -2^P g v g (22) 

§- t (P S E) = ~[(37-l)p g 6 + 2p gV g (23) 

^(p g e) = -^(3 7 -l)p g e. (24) 

We solve the radiative cooling term using a semi-implicit numerical ODE integrator. For stability 
both operators require timestep limitations: for the cosmological terms, we limit the change in a(t) 
during the timestep to 10% or less, while for the radiative cooling term we use 

At cool = 0.25^ . (25) 

These timestep restrictions are considered along with the hydrodynamical and iV-body restrictions 
in determining the actual timestep via equation (20). 

We also treat the gravitational acceleration terms in an operator-split fashion, but because 
they involve spatial derivatives, we must make a small modification to the PPM method to ac- 
comodate them. This modification is described by Colella & Woodward (1984). In Section 3.2.1 
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we discuss further the discretization of the Poisson equation and the expressions used for the cell- 
averaged gravitational acceleration in PPM. We impose the following timestep constraint due to 
the gravitational acceleration gijk- 



At 



grav 



mm 



Ax,; 



Ay,- Az k 



,1/2 



§ rid LlSaw'fcl ' \dy,ijk\ 
In practice this criterion rarely determines the timestep. 



\dz,ijk\_ 



(26) 



3.1.2. The piecewise-parabolic method 

To calculate the time-averaged fluxes we use the piecewise-parabolic method (PPM; Colella 
& Woodward (1984)), a high-order Eulerian generalization of Godunov's (1959) scheme. Godunov 
schemes achieve very good resolution of flow discontinuities such as shocks without excessive dis- 
sipation or oscillations by solving a nonlinear flow problem at each cell interface. This property 
makes PPM ideal for solving a variety of astrophysical flow problems, since such problems often 
require narrow hydrodynamical features to be resolved with a small number of cells. 

Since the PPM algorithm is discussed in detail elsewhere, we limit ourselves to describing the 
special features of the version of PPM used in COSMOS. We implement PPM using the direct 
Eulerian formulation, and we include the modifications developed by Colella Sz Glaz (1985) for 
handling general equations of state using a variable ratio of specific heats 7. We use the approximate 
Riemann solver described by van Leer (1979); rarefactions are locally approximated as shocks, 
following Colella (1985). 

PPM uses information from four cells on either side of each cell at each timestep. At the 
boundaries of the computational region we therefore maintain four 'ghost zones' containing bound- 
ary information. We set the values in these zones before applying each ID differential operator. 
Periodic boundaries, which are appropriate for large-scale structure simulations, are implemented 
by setting the ghost zones equal to the corresponding interior zones on the opposite side of the grid. 
For galaxy cluster evolution problems it is more desirable to control inflow and outflow rates than to 
have matter which disappears from the grid reappear on the opposite side. For inflow boundaries 
one sets the ghost zones to prescribed values; for outflow, one typically uses Neumann or zero- 
gradient boundaries (Roache 1998). In the latter case the interior zones adjacent to the boundaries 
are simply copied into the ghost zones. For most astrophysical outflow applications this approach is 
fine, although strictly speaking when the outflow is subsonic one should use a method-of-lines-based 
approach to prevent reflections from propagating into the interior (Thompson 1987). 

However, when zero-gradient outflow boundaries are used with a self-gravitating fluid, one 
must take care that the flow in the interior zones adjacent to the boundary is indeed directed 
outward, else the copying of the velocity component normal to the boundary into the ghost zones 
will create an artificial (and destabilizing) inward flux. For isolated problems we do not prescribe 
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an inflow pattern. In such cases we therefore implement a form of vacuum boundary in which 
Neumann conditions are used unless the normal velocity in the zones just interior to a boundary 
is directed inward. If the normal component of velocity is directed inward, we use Dirichlet (zero- 
value) boundaries for it and Neumann conditions for the other variables, and we zero the external 
fluxes (resulting from the solution of Riemann problems on the external boundaries) if they happen 
to be inwardly directed. In isolated calculations, therefore, matter which leaves our computational 
grid disappears altogether and cannot return. Since the gas is self-gravitating, this would cause 
problems with the potential if a large clump of material were to leave the grid suddenly. We must 
therefore take care to select a large enough grid to prevent substantial amounts of material from 
leaving during the course of a simulation. 



3.2. Gravitation 

3.2.1. Finite-volume discretization 

At the end of each timestep we solve the Poisson equation (5) for the combined gravitational 
potential of the gas and dark matter. If we apply the spatial convolution operator to this equation 
and use Taylor expansions about the cell walls to obtain cell-wall averages of d(p/dx, d(f)/dy, and 
d(p/dz, we obtain the following second-order discretized version, in which only the x-derivatives are 
presented for clarity: 

B £,i<t>i+2,jk + B ti<f>i+i,jk + {B^i - B'^ijk - B'^i-xjk - B'^i-zjk + ...= (27) 

4vrG 



o 3 



[{Pg,ijk + Pdm,ijk) - Pg + PdmJ 



Here we define 



At = ±- 



' Axi(Axi + Ax i± i + Axi± 2 ) 
? ± _ A ±2 Axj±! + Ax i± 2 
Axi + Ax i± i 



B£i - "A a a ( 2 8) 



M ~ i {Axi + Ax i± i)(Ax i± i + Axi± 2 ) 
B ± = A ± Axi-Ax i±1 



Ax i± i + Axi± 2 



The spatially averaged total density (p g + pdm) is subtracted from the total density to make the 
potential tend to zero at large distances from a point source. The cell-averaged dark matter density 
Pdm,ijk is provided by the particle-mesh code. The complicated form of equations (27)-(28) is due 
to the nonuniform grid. If the spacings Axi are identical, the left-hand side of equation (27) reduces 
to the usual second-order difference (<^+i — 2(f>i + 0j_i)/Ax 2 . 

The densities on the right-hand side of equation (27) are cell-averaged quantities, so the po- 
tential which results from solving this equation is also a cell-averaged quantity. This raises the 
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question of how to obtain the gravitational acceleration g^, which the discretized gas momentum 
and energy equations require to be a cell-averaged quantity since it is a source term. With the 
spatial convolution operator we can obtain an expression for the x-component of gijk'- 



9x,ijk 



f-Xi+Axi/2 r-yj+Ayj/2 f-z k +Az k /2 o 

/ / / dV-f^x') (29) 

Jxi-Axi/2 Jyi-Ayj2 J z k -Az k l2 ox 



Axj Ayj Az k J Xi -Axi/2 J yj -A yj /2 Jz k -Az k /2 
= [^+1/2 ,3k ~ <t>i-l/2,jk] ■ 

We compute the cell-wall-averaged potential {<Pi+i/2jki etc.) from the cell averages <pijk using Taylor 
expansions of 0(x) about the cell walls in each direction. As an example, the results for 4>i-i/ 2 jk 
and <p i+ i/ 2 ,jk are 

0i±i/2jfe = C^i4>ijk + Cf ti 4>i±ijk + C^i4>i±2,jk , (30) 



where 



n ± _ Azy-i(AcCi±i + Axi± 2 ) 



y 0,i 



(Axi + Ax i± i)(Axi + Ax i± i + Axi± 2 ) 



c ± = Axi [Ax i ± 1 {2Axj + 3(Aa i± i + Ax i±2 )) + Axi± 2 {Axi + Aa i±2 )] , . 

M ~ (Axi + Ax i± i)(Ax i± i + Ax i±2 )(A Xi + Ax i±1 + Ax i±2 ) 
n ± _ AxjAxj-ti 



/2,i " (Aar a . ± i + Ax i±2 )(Ax i + Ax i± i + Ax i±2 ) ' 

For each cell we calculate gjj^ in this way using the cell average of cf> in that cell and its four nearest 
neighbors in each coordinate direction, yielding second-order accuracy for gijk- If the mesh spacing 
is uniform, the resulting expression for g x ,ijk is 

9x,ijk = [ ( t ) i+2,jk - 5(4>i+l,jk - 4>i-l,jk) - 4>i-2,jk] ■ (32) 



3.2.2. The multigrid solver 

We solve the Poisson equation (5) using the full multigrid (FMG) method (Brandt 1977; Hack- 
busch 1985). This method is as fast as the direct transform-based methods but can be parallelized 
more easily and in a more machine-independent fashion. It also generalizes more easily to nonuni- 
form grids and nonperiodic boundaries. In addition to the finite-volume discretization described 
in the previous section, our multigrid implementation features the capability of handling isolated 
boundary conditions. Note that our use of multigrid to compute the gravitational potential is to 
be distinguished from the solution of potential-flow problems in large-scale structure simulations, 
in which the hydrodynamical equations also are solved using multigrid techniques. 

To implement FMG, we begin by constructing a hierarchy of nested grids, each of which is 
twice as coarse as the previous one (Figure 1). The hierarchy starts with a grid only a few zones 
across, on which the Poisson equation can be solved directly, and ends with a grid identical to 
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that used for the hydrodynamical equations. On each grid level FMG applies an iterative solution 
method (here, Jacobi), bringing errors on all length scales into convergence at the same rate. 

Our multigrid scheme uses an ascending V-cycle (Figure 2). We begin with a guess for the 
solution on the coarsest grid. We first apply a few Jacobi iterations to the guess (this smooths the 
error modes with wavelengths which are short relative to the grid spacing), then 'prolongate' this 
approximate solution to the next finer level. Our prolongation operator is simple trilinear inter- 
polation using cell-averaged quantities with nonuniform (though nested) grids. After performing a 
few Jacobi iterations at this level, we then calculate the 'residual', the difference between the guess 
and the result of applying the Laplacian operator to the density on this grid. Since the Poisson 
equation is linear, the solution to the residual equation is just the correction which must be added 
to the original guess to obtain the true solution. However, since long-wavelength error modes are 
responsible for the Jacobi method's slow convergence, the residual on the second level will be dom- 
inated by errors that have long wavelengths relative to this grid and shorter wavelengths relative 
to the coarsest grid. To remove these errors, then, we 'restrict' the residual back to the coarsest 
grid and there solve exactly for the correction (actually several Jacobi iterations on this level are 
sufficient). Since we are using finite volumes, restriction is equivalent to simple averaging over the 
fine-grid zones which lie within a given coarse-grid zone. After solving on the coarsest grid for the 
correction, we prolongate the correction back to the second level and apply it to the solution guess 
on this level. We repeat this V-cycle twice, then prolongate the solution guess to the next finer 
level, where we perform three more V-cycles, each time restricting all the way back to the coarsest 
level. We repeat this process until we reach the finest level (the original grid). 

Boundary zones for the multigrid solver must be maintained both on the sides and at the 
corners of the computational cube for each level of the grid hierarchy, although the corners are 
only used for prolongations. The grid levels are nested, so their external boundaries coincide, and 
Dirichlet, Neumann, or periodic boundaries are simple to implement. The primary complication 
lies in distinguishing between cell-averaged values and cell-wall-averaged values of the potential 
and the density; the cell averages in the ghost zones must be chosen so as to make the cell-wall 
averages of the potential on the boundary possess the desired properties. We do this by writing 
the potential as a Taylor series about the boundary location in the direction perpendicular to the 
boundary, then average the series over nearby cells, obtaining a set of equations for the cell-wall 
averages of the potential and its derivatives in terms of the (known and unknown) interior and 
ghost-zone cell-averaged potential. We invert these equations (up to some limiting order) to obtain 
the unknown ghost-zone cell average(s) in terms of the known boundary values and interior-zone 
cell averages. 

Isolated boundaries are boundaries outside of which the source function is identically zero. 
Implementing them requires that we specify the value of the potential on the boundary surface. In 
order to solve for the gravitational potential of isolated matter, we use a variant of James' (1977) 
method. We first calculate the potential assuming Dirichlet boundaries. We then find the image 
distribution required to make the potential identically zero outside the boundary by evaluating the 
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Laplacian of this Dirichlet solution on the boundary. Finally, we calculate the isolated potential 
of the (hollow) image distribution by computing its moments up to some maximum multipole 
and solving with boundary values appropriate to the corresponding multipole expansion of the 
potential. By subtracting this from the Dirichlet solution we obtain the desired isolated potential. 
We reduce the number of multipole moments required by first finding the center of mass of the 
image distribution and then performing the multipole expansion about this point. Our method is 
somewhat similar to that of Miiller & Steinmetz (1995), who use a spherical-harmonic expansion 
of the original source function about the center of a spherically symmetric grid to compute the 
potential at each of the interior grid points. 



3.3. Dark matter 

To handle collisionless components, such as dark matter and stars, we use an iV-body code 
based on the particle-mesh (PM) method (Efstathiou et al. 1985; Hockney & Eastwood 1988). This 
technique speeds force calculations for particles by converting particle positions to densities on a 
grid, then solving the Poisson equation on the grid using a fast direct or multigrid solver, and 
finally interpolating the potential from the grid to obtain forces at the particle locations. Since we 
already need to solve the Poisson equation on a grid for our hydrodynamical code, the particle-mesh 
method is ideally suited for integration with this code. We simply add the equivalent densities for 
the particles to the grid densities for the gas, then solve the Poisson equation as usual. In our case, 
grid densities for the particles are computed using the cloud-in-cell (CIC) operator. We also use 
this operator to compute interpolated forces. 

Other iV-body methods, including particle-particle-particle-mesh (P 3 M) and the tree methods, 
are often used to follow the dark matter in cluster simulations. The P 3 M method (Efstathiou & 
Eastwood 1981) extends the resolution of particle-mesh by using direct summation for the force 
between particles lying within a single zone. Tree methods (Hernquist 1987), on the other hand, 
approximate the force due to distant groups of particles using their low-order multipole moments. 
However, because we are using a grid-based hydrodynamical method, we do not expect to resolve 
features smaller than a single zone. Therefore the considerable execution time expended in the 
direct summation part (which scales as the square of the number of particles) of P 3 M would be 
wasted. The extra resolution provided by tree codes would likewise be wasted. Particle-mesh has 
another advantage for us over tree methods in the ease with which it can be integrated with the 
gas code. 

The main difference between the PM part of COSMOS and the 'classic' PM codes (e. g., 
Efstathiou et al. 1985) is our use of a variable timestep. In order to keep the computation accurate 
to second order, we need a more complicated version of the standard leapfrog update of position 
and velocity which makes use of the gravitational acceleration on each particle stored from the 
preceding timestep. The resulting update steps (given here only for the x component of position 
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and velocity) are 
for the position and 
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for the velocity. Here we define A = — 2a/ a. For the first half-timestep we use a first-order accurate 
expression to obtain vi/ 2 from the initial conditions. Note that the derivatives of the potential from 
the previous timestep (here denoted V</> n _i) must be retained. This is also necessary in the PPM 
code in order to make second-order accurate gravitational corrections to the Riemann solver, as 
described by Bryan et al. (1995). 



4. Code tests 

As with any scientific code of this complexity it is necessary to apply COSMOS to a number 
of test problems with known solutions before using it on research problems. This is important not 
simply for the purpose of verifying that the code works as it should, but also because it gives us an 
intuitive understanding of the code's strengths and weaknesses. Such an understanding is essential 
for making effective use of any hydrodynamical code. Accordingly, we have performed a number 
of tests which exercise the various code modules in different combinations. We describe the test 
problems and the code's performance on each in this section. 



4.1. Hydrodynamics tests 

4-1.1. Sod shock-tube problem 

The widely used test problem due to Sod (1978) is a Riemann shock tube problem with 
initial conditions specially chosen to produce all three types of fluid discontinuity (shock, contact 
discontinuity, and rarefaction). In this test we create two regions of constant density and pressure, 
separated by a plane whose normal forms specified angles with the x and y axes. To the left of 
the plane we set pl = PL = 1, and to the right we set pr = 0.125 and pr = 0.1. The ratio of 
specific heats, 7, is 1.4, and the fluid is everywhere at rest. This test is purely hydrodynamical; 
no gravitational potential is used. The Sod test enables us to determine if our code satisfies the 
shock jump conditions and whether it correctly determines the speed of each nonlinear wave. By 
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performing this essentially one-dimensional test in three dimensions at an angle to the grid, we can 
also determine how well we can resolve shocks for realistic grid sizes in 3D. 

Figure 3 plots our numerical solution to the Sod problem against the analytic solution at 
t = 0.206. We used a 64 3 uniform grid in the box [0, l] 3 ; the figure shows profiles taken as a 
function of distance s along the line segment connecting (0,0,0) and (1,1,1) (normal to the shock 
plane). The initial fluid discontinuity was located at s = 0.577. In constructing the initial conditions 
we smoothed the discontinuity in each variable over one zone's width to reduce the starting error 
resulting from the presence of all three nonlinear waves within one zone at early times. Since errors 
in the contact discontinuity propagate at the speed of the fluid (and hence of the discontinuity), 
they accumulate at the discontinuity, and any starting error will therefore be present at all later 
times. Even with the initial smoothing a 2% error is present in the density in the three zones 
in front of the contact discontinuity; this results in a similar error in the specific internal energy. 
Apart from this error and a slight underestimation of the slopes in the leftward-moving rarefaction 
wave, the numerical solution is nearly exact. The position of each of the three discontinuities is 
correct, and the shock and contact discontinuity are each resolved using about two zones, despite 
the fact that the shock is moving at an angle to the grid. This excellent shock handling is one of 
the primary reasons for using PPM. 



4-1-2. Sedov explosion problem 

To determine how well our hydro dynamical code respects the 90° rotational symmetry of our 
grid, and to determine whether this behavior carries over to a nonuniform grid, we studied the 
expansion of a strong spherical shock wave into a uniform medium. In this problem, which is again 
purely hydrodynamical, we deposit a quantity of energy E = 1 into a small sphere of radius Sr at 
the center of the grid. The pressure inside this volume, p' , is given by 

where for this problem we use 7 = 7/5. Everywhere the density is set equal to po = 1, and 
everywhere but the center of the grid the pressure is set to a small value, po = 10~ 5 . The fluid 
is initially at rest. Sedov (1959) first obtained a self-similar analytic solution to this problem. 
A spherically symmetric shock wave develops; the density, pressure, and radial velocity are all 
functions of 

t = m = c(w) ' (36) 

where C is a numerical constant depending only on 7. Just behind the shock front at £ = 1 we 
have 
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where u = dR/dt is the speed of the shock wave 
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In Figure 4 we plot the analytic solution at t = 0.0508 against the angle average of the numerical 
solution found by COSMOS using a uniform 64 3 grid in the box [0, l] 3 . At this time the shock front 
is located at r = 0.314. Since the grid is Cartesian, in the initial conditions we have attempted to 
minimize geometrical effects due to the shape of the grid by using an initial sphere of radius (8r) 
3.5 zones. For each zone containing part of this sphere we weight the initial pressure according to 
the fraction of the zone which lies within the sphere using Monte Carlo averaging. Despite these 
efforts some small errors of geometrical origin are still present, particularly in the velocity field, 
where some oscillation (~ 7%) can be seen behind the shock. The position of the shock itself is 
accurate to 1-2 zones, which is also the width of the shock in the numerical solution. However, 
because the shock is narrower than one zone in the analytic solution, the density and pressure 
in the numerial solution do not reach their maximum analytic values p\ and p±. The density is 
also slightly underestimated near r = 0.25, as is the central pressure. Given the use of a uniform, 
relatively coarse Cartesian grid for this spherically symmetric problem, the position and shape of 
the shock are well-determined. 

To determine the effect of a nonuniform grid on our Sedov solution, we solved the same problem 
on grids with different degrees of nonuniformity. (The grid in question is static, not an adaptive 
grid designed to track the explosion more accurately.) In these runs the x-axis is nonuniform, and 
the y and z axes are uniform. Along the x-axis the innermost -/V; nner zones are uniform with spacing 
A. Outside this region the zones increase in width by rj per zone toward the outer edges. Each 
axis uses 64 zones; the y and z axes each enclose the range [0, 1], and for the x-axis we fix iVj nner 
at 16 and A at 1/64, allowing rj to take on several values between and 20%. (Hence the box 
size in x also varies.) For each run the explosion is thus allowed to develop up to a certain point 
within the same uniformly gridded region; thereafter part of the shock front expands into a region 
of nonuniform gridding. 

Figure 5 shows, as a function of radius, the percentage change in the angle-averaged density at 
t = 0.0508 for several values of r\ up to 20%. Each density profile has been interpolated to the same 
uniform grid. Because the shock is only partly resolved even on the uniform grid at t = 0.0508, 
the increased zone width for large r\ forces the shock to be substantially broadened, leading to 
an overestimate of the density ahead of the shock and an underestimate behind it. The average 
density is also underestimated by a constant 1-2% well behind the shock for rj > 10% because the 
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broadened shock does not sufficiently compress the ambient medium. For 77 = 5% the magnitude 
of these effects is much smaller. 

As the explosion proceeds in the analytic solution, the width of the shock front increases, and 
its amplitude decreases. If the zone size increases much faster than the rate of increase of the width 
of the shock, the numerical solution will continue to degrade as the shock propagates outward. 
For this problem the value of r] at which the zone size begins to increase faster than the shock 
width appears to lie between 5% and 10%. In general, the amount of nonuniformity to use must 
be balanced against the need to resolve fine flow structures near the edge of the grid. If the zone 
size is permitted to be too large, the zone-averaged density, pressure, and velocity may be correct, 
but the position and size of flow structures will be poorly determined, and their dynamical effects 
(such as heating and compression) on the surrounding gas will also be in error. 



4-1.3. Jeans instability problem 



Using the Jeans instability we verified that the gas dynamical terms involving the gravitational 
potential are correctly implemented in COSMOS. To do this we studied the dispersion relation of 
stable perturbations to a uniform, self-gravitating medium. The initial conditions for this problem, 
which we solved in two dimensions with periodic boundaries, are, at time t = ti, 
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(39) 



where the perturbation amplitude S <C 1, and where oj is defined as follows. The stability of the 
perturbation is determined by the relationship between the wavenumber k = \[2k and the Jeans 
wavenumber kj (Chandrasekhar 1961), where kj is given by 



and where cq is the sound speed: 
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If k > kj, the perturbation is stable and oscillates with frequency 

uj = \Jc 2 Q k 2 -4irGp ; 
otherwise it grows exponentially, with a characteristic timescale given by r = (iuj)^ 1 . 
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We checked the dispersion relation (42) for stable perturbations with 7 = 5/3 by fixing po and 
Po and performing several runs with different k. We followed each case for roughly ten oscillation 
periods using a uniform 64 2 grid in the box [0, y / 7T7/2] 2 with units in which po = po = G = 1. 
(The box size is chosen so that kj is smaller than the smallest nonzero wavenumber which can 
be resolved on the grid. This prevents numerical errors at wavenumbers less than k j from being 
amplified by the physical Jeans instability.) We then computed the oscillation frequency in each 
case by measuring the time interval required for the density at the center of the grid to undergo 
exactly ten oscillations. The resulting dispersion relation is compared to equation (42) in Figure 
6. It can be seen from this plot that our code correctly reproduces equation (42). At the highest 
wavenumber (k = 62.132), each oscillation is resolved using only about 9 cells, and the average 
timestep (which depends on cq, Ax, and Ay, and has nothing to do with k) turns out to be about 
one-fifth of an oscillation. Hence the frequency determined from the numerical solution for this 
value of k is somewhat more poorly determined than for the other runs. Overall, however, the 
frequencies are correct to about 1%. 



4- 1-4- Z el' dovich pancake 

The cosmological pancake problem (Zel'dovich 1970) provides a good simultaneous test of the 
self-gravity and cosmological expansion code. Analytic solutions well into the nonlinear regime 
are available for both iV-body and hydro dynamical codes (Anninos & Norman 1994), permitting 
an assessment of the code's accuracy. After caustic formation the problem provides a stringent 
test of the code's ability to track thin, poorly resolved features and strong shocks using most 
of the basic physics needed for cosmological problems. Also, as pancakes represent single-mode 
density perturbations, coding this test problem is useful as a basis for creating more complex initial 
conditions. 

We set the initial conditions for the pancake problem in the linear regime using the analytic 
solution given by Anninos and Norman (1994). In a universe with Qq = 1 at redshift z, a pertur- 
bation of wavenumber k which collapses to a caustic at redshift z c < z has comoving density and 
velocity given by 
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where p is the comoving mean density. Here x e is the distance of a point from the pancake midplane, 
and X£ is the corresponding Lagrangian coordinate, found by iteratively solving 

l + z c smkxe 

x e = xe — ; . (44) 

1+z k y ' 

The temperature solution is determined from the density under the assumption that the gas is 
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adiabatic with ratio of specific heats 7: 

T(x e ; z) = (l + z) 2 f M 



1 + ZficA p(x e ;zM) 



l + z 
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(45) 



The mean temperature Tg^ is specified at a redshift z^. 



At caustic formation (z = z c ), planar shock waves form on either side of the pancake midplane 
and begin to propagate outward. A small region at the midplane is left unshocked. Immediately 
behind the shocks, the comoving density and temperature vary approximately as 
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At the midplane, which undergoes adiabatic compression, the comoving density and temperature 
are approximately 
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To test the convergence of the code, we performed several different one-dimensional pancake 
test runs with fixed perturbation wavelength and caustic redshift and varying spatial resolution. 
Each run used k = 27r/(10 Mpc), z c = 5, and an initial redshift Zi = 50. In each run we also fixed 
z M = 200, T fid = 550 K, 7 = 5/3, ft = 1, and H = 50 km s" 1 Mpc" 1 . We assumed a 75% 
H/25% He composition for the purpose of computing the pressure using the temperature. Using 
the analytic solution (eqs. [43] - [45]), we computed LI error norms for density and temperature 
in each run at z ~ 7, in the mildly nonlinear phase of the collapse. We define the LI error norm 
as in equation (50) below. The results appear in Figure 7. The density converges slowly, with 
the error varying approximately as Ax ' 6 . The temperature converges more rapidly, with the error 
varying as Ax 3 ' 5 for Ax > A/128 and then as Ax 1 ' 3 for smaller Ax. The roughly linear asymptotic 
convergence for the temperature is consistent with other work, for example that of Bryan et al. 
(1995), but the density convergence rate is somewhat slower. Nevertheless, the absolute error of a 
few percent for Ax ~ A/100 — A/200 is consistent with their results for both variables. 

To examine the performance of COSMOS on this problem in multidimensions, we performed 
a run with the same initial conditions using a doubly periodic 256 2 grid. The pancake midplane 
was inclined 45° to the x axis, and the zone spacing was chosen so that two wavelengths fit into the 
box diagonal. In Figure 8 the results of this run are compared to the 256-zone ID solution at three 
different redshifts corresponding to the mildly nonlinear regime [z = 7.4), the time immediately 
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following caustic formation (z = 4.4), and a time well into the nonlinear, post-caustic evolution 
(z = 0). The profiles for the 256 2 run are derived by interpolating along the grid diagonal. In this 
run only 128 zones fit within a perturbation wavelength, so the spatial resolution is one-half that of 
the ID 256-zone run. With the exception of the innermost 1-2 zones, the 2D results are typically 
within a few percent of the ID results. Differences of note include a 20% underdensity in the 2D 
case at z = 7.4 and a significant overestimate in the central temperature at z = 4.4. In the latter 
case the caustic is unresolved by the mesh, and only the innermost zone is hot. The velocity in the 
innermost zone is different by as much as a factor of two at all three redshifts. 



4.2. Tests of the Poisson solver 



4-2.1. Miyamoto-Nagai potential 



We tested the potential solver with isolated boundaries using the Miyamoto-Nagai potential 
(Miyamoto & Nagai 1975). This is an axisymmetric, flattened potential designed to mimic the light 
distribution of a disk galaxy. It is given in cylindrical coordinates by 
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where M is the total mass, and the ratio of the axis parameters b/a determines how flattened 
the potential is. As b/a — ► 0, 0mn tends toward the Kuzmin disk potential, and as b/a — ► oo, 
4>m~n tends toward the spherically symmetric Plummer potential. Thus for different b/a values </>mn 
contains different contributions from high-order multipole moments, making this a good test of our 
moment-based isolated boundary solver. The density function corresponding to this potential is 
straightforward to derive using the Poisson equation; it is 
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Using this density function we computed Miyamoto-Nagai potentials with the isolated Poisson 
solver in COSMOS for b/a = 0.1, 1, and 10 using a 128 2 x 64 grid in boxes containing 98-99% 
of the total mass. The innermost one-half of the zones in each dimension were uniform, and the 
remainder increased in width by 5% per zone toward the outside edges. We used units in which 
G = M = a = 1. 

In our treatment of isolated boundaries we compute boundary values for the image potential 
using a truncated multipole expansion of the image mass distribution. To verify that this expan- 
sion was implemented correctly, we performed runs with four different values of £ max , the largest 
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multipole moment. In Figure 9 we plot average errors for our computed potentials as a function of 
4nax for our three values of b/a. We define the error E(£ max , b/a) as 



where the sum is taken over all of the cells in the grid, and the total volume V of the grid is 
the sum of cell volumes AV^- Figure 9 shows that including only the first three multipole terms 
already gives a fairly good estimate of the potential; the maximum percentage error in any cell for 
^max = 2 is between 6.5% and 7.5%, and the plotted averages range from 2.2% to 2.7%. Increasing 
the value of £ m &x brings the average and maximum errors down substantially, showing that the 
isolated boundary solver is including higher-order moments properly. The curves begin to level off 
at values less than 1% above £ max = 10; this level of error is consistent with the amount of mass 
which lay outside the grid and hence was excluded from each calculation. The maximum error at 
4nax = 14 is about 1.5% in all three cases. 

For other problems the multipole content of the density field may differ, so these test results 
do not show that 1% errors will be achieved with. f max f» 10 for any arbitrarily chosen density field. 
However, this will be true even for quite flattened density fields. Errors caused by underestimating 
^max will be largest near the boundaries of the grid as long as most of the mass is near the center. 



Because of the particle-mesh force-smoothing procedure, dark matter particles experience an 
effective potential that deviates from the Newtonian 1/r dependence at small interparticle sepa- 
rations r. This deviation is greatest on length scales less than or comparable to the zone spacing 
Ar. In addition, the introduction of a grid breaks the rotational symmetry of the equations of 
motion, so interparticle forces are not isotropic for r <^ Ar. When performing simulations in which 
small-scale structure (relative to the grid in use) appears, it is important to have an estimate of 
the magnitude of these errors and the value of r/Ar at which they become important. 

To quantify these effects, we performed a test in which we chose 10,000 randomly placed pairs 
of particles on a 32 3 grid, computed forces for each pair using the multigrid solver and particle-mesh 
code from COSMOS, and tabulated the resulting forces as functions of the particle positions. The 
first particle in each pair was chosen to lie within one zone of the center of the grid, while the second 
particle was chosen to lie between 0.1 and 10 zones from the first, with a random position angle 
relative to the first. In this test each particle had unit mass, and the value of G was unity. For 
the ith particle of the j'th pair we computed radial (/ij, r ) and azimuthal {fijfi) force components 
relative to the line connecting the particle to its partner; then for the pair we computed the average 
of the magnitudes of the radial and the azimuthal components for the two particles: 




(50) 



4-2.2. Particle-mesh force resolution test 



fj,x = (\fij,x\ + \f2j,x\) 



(51) 



- 22 - 

where X = r, 9. We collected the particle pairs in radial bins of logarithmic width 0.05, then 
computed mean and RMS deviations for the binned forces. 

In Figure 10 we show the results as functions of radius. Overall the results are comparable to 
those obtained by Efstathiou et al. (1985): our effective particle force resolution is about two zones. 
For r = Ar the mean radial force is 50% of the expected value, and the force anisotropy (ratio of 
mean azimuthal to mean radial force) is about 13%. For particle separations smaller than this, the 
mean radial force is proportional to r, and the force anisotropy is roughly constant at around 7%. 
As r — > oo the azimuthal force drops as r~ 3 , so the force anisotropy drops approximately as r -1 . 
At r f» 9 the azimuthal force turns upward slightly; this illustrates the effect of placing the second 
particle of some pairs within 4-5 zones of the grid boundary. We discuss this effect in more detail 
below. 

For particle separations less than one zone, the RMS deviations in the radial and azimuthal 
force components are each roughly constant at 20% and 5%, respectively, of the mean radial force. 
For r > 2Ar the radial RMS deviation drops to about 7% of the radial mean, while the azimuthal 
RMS deviation decreases almost as fast as the mean azimuthal force, becoming less than 0.2% of 
the radial mean at r = lOAr. The azimuthal RMS deviation also shows a slight upturn for particles 
close to the boundary. 

The deviations from a 1/r 2 force law at large r result from a self-force felt only by particles 
close to the edge of the grid. As Hockney and Eastwood (1988) show, self-forces arise in PM 
schemes when either the density assignment and force interpolation operators differ or the effective 
Green's function Gijkimn fails to satisfy the criterion 

Gijklmn = ^Imnijk 

The effective Green's function is related to the zone-averaged gravitational field by 

= ^ijklmnP^mn^Vlmn j (53) 
Imn 

where pfjj^ is the density assigned to the grid by the particle-mesh scheme, AV/ mn is the volume 
of a single grid zone, and the sum is taken over all zones. Because we use the CIC operator both 
to assign densities to the grid and to interpolate forces to particle positions, any self-forces are due 
to errors in the potential which cause G to violate equation (52). In our case, errors in G arise 
because we use a finite multipole expansion to obtain boundary conditions for isolated problems. 

We have examined this effect by performing another test in which one particle is placed at the 
center of the grid and the second particle is placed at increasing distances from the center along 
one coordinate axis. By varying the grid size, we can use this test to determine how close to the 
edge of the grid a particle can be before it experiences a significant self-force. For this test we use 
4nax = 10 in the multigrid solver. Figure 11 compares the results on grids with 16 3 , 32 3 , and 64 3 
zones. In each case the force on the second particle begins to deviate from the expected value when 
the second particle is about five zones from the edge of the grid. 



(52) 
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Note that the two-particle potential is more susceptible to this effect than the Miyamoto-Nagai 
test potential because of the relative importance in this case of mass close to the edge of the grid. 
The effect is important only for small numbers of particles close to the edge of the grid; it is not 
important for particles near the boundary orbiting in the potential of a large group of particles 
near the center, and it is not present when periodic boundary conditions are used. 



4.3. iV-body tests 

4-3.1. Time integration accuracy 

To account for the varying timestep size, COSMOS implements equation (34) to update the 
dark matter positions and velocities. To test this algorithm, we set up two particles in a simple 
circular orbit and monitored the accuracy of the orbit as we varied the timestep At and the zone 
size Ax. In addition to testing the algorithm, this problem enables us to determine the regions in 
(Ax, At) space that are appropriate for larger, more interesting problems. 

Figure 12 shows the root mean square (RMS) error in the interparticle separation after six 
orbital periods as a function of Ax. For Ax larger than the particle separation, the 1/r 2 force is 
calculated incorrectly due to the mesh smoothing. The problem of resolving forces on scales smaller 
than a zone size is a well-known one; Figure 12 simply points out that our force resolution is roughly 
one or two zones. For very small Ax, on the other hand, each particle traverses more than one 
zone in a single timestep, which again leads to large errors. To alleviate this second problem, in 
production runs we choose At with the requirement that no dark matter particle travel more than 
a fraction a^m of a zone during a single timestep: 



Atdm = 0dm min 

particles 



Axj Ayj Az k 



p 1 V ' V 
Va Va V a 

am,x dm, y dm, z 



(54) 



where (i,j, k) are the indices of the zone containing the pth particle. Figure 12 shows that setting 
this fraction equal to one is an acceptable choice. The two groups of points in this plot correspond 
to different ratios of the timestep to the orbital period. If we denote this ratio by /, then we require 

for accuracy. The squares indicate runs with / = 0.046, and the triangles indicate runs with 
/ = 0.015. If we choose a^m = 1) then the two groups of runs require Ax/2r > 0.14 and 0.047 
respectively. For the first group, runs with Ax/2r > 0.1 are observed to be relatively accurate, 
while for the second group, the error is low for Ax/2r > 0.06. In the regime where the zone size 
is not too big or too small relative to the timestep, we see that the cumulative error does indeed 
scale as At 2 : the algorithm of equation (34) is accurate to second order. 
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4-3.2. Spherical collapse problem 



A well-known cosmological problem with an analytic solution is the case of a spherical over- 
density (e. g., Padmanabhan 1993). The solution to this problem is best described in terms of the 
radius R(t) which encloses a mass M at time t: 



R(t; M) = R(U;M) sin 2 



l + Si 



(56) 



Here U and Si denote the time and fractional overdensity at which the simulation starts. The 
parameter 8 is a timelike parameter; in a flat, matter-dominated universe, the redshift can be 
expressed in terms of as 



1 + Z = (1 + Zi) 



1 + 7^4^ ( 0-81110-0* + sinfl, 

4 ^3/2 



-2/3 



(57) 



Equation (56) holds until the radius reaches a maximum value at = w, after which time virializa- 
tion occurs. 

On a 32 3 grid with 32 3 particles we set up spherically symmetric initial conditions. In units in 
which the cell size is unity, within a radius of five units from the center of the grid we constructed 
a uniform overdensity with Si = 0.3. In an annulus ranging from r = 5 to r = 8 we placed 
an underdensity so that the average density within r = 8 was equal to the background density. 
Figure 13 shows the time evolution of several (comoving) radii, compared with the analytic result 
of equation (56). The larger shells track the analytic solution well, verifying that all factors of a(t) 
are correctly implemented in the dark matter code. As expected, the small shells do worse: the 
force on these shells comes predominantly from particles concentrated at a distance of one to two 
zones. The deviations from equation (56) for these shells are consistent with the underestimate of 
the gravitational force introduced by the particle-mesh smoothing at this distance. 



4-3.3. One-dimensional Plane Wave 

The Zel'dovich pancake (Section 4.1.4) also serves as a test of the iV-body code. To initialize 
this problem, we place particles at the zone centers of a 16 3 -zone grid, then perturb their positions 
slightly in the x-direction by amounts 

Sxi oc cos(2irkxi/L) (58) 

with Xi the unperturbed positions and k an integer. Small-scale perturbations correspond to large 
k. We consider only perturbations which vary on scales larger than a zone width, ie. those with 
k « N x . If the particles' peculiar velocities Sv x ^ are also proportional to cos{2irkxi / L) , then in 
the linear regime (assuming a flat universe) both Sxi and Sv Xj i should grow with time as the scale 
factor a(t). 
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To determine how accurately their code could follow such perturbations, Efstathiou et al. (1985) 
tabulated the RMS deviation of position and velocity from the analytic solution for particles on 
a 16 3 mesh. Since we solve Poisson's equation differently than the standard PM code, it is useful 
to compare our results with theirs. Figure 14 shows the RMS errors in position and velocity for 
linear perturbations with several different values of the wavenumber k. The multigrid solver and 
variable-timestep integrator used in COSMOS produce results quite similar to (but slightly better 
than) the FFT/fixed-timestep technique used by Efstathiou et al. (1985). 

5. Conclusion 

In this paper we have described COSMOS, a new hybrid code for solving cosmological problems 
involving gas and collisionless matter, including self-gravity, cosmological expansion, and radiative 
cooling. COSMOS solves the inviscid fluid equations using the piecewise-parabolic method on 
a static, nonuniform structured grid. The code treats collisionless matter using the cloud-in-cell 
variant of the particle-mesh method, and it computes the gravitational potential using a linear full 
multigrid solver. COSMOS supports both periodic boundaries, suitable for large-scale structure 
problems, and isolated boundaries, suitable for simulations of isolated systems. This code has 
already been used to study mergers between clusters of galaxies by Ricker (1998) and Ricker & 
Sarazin (2000). 

Several new features of COSMOS distinguish it from other existing PM-PPM codes. In par- 
ticular, the multigrid Poisson solver can determine the gravitational potential of isolated matter 
distributions and properly takes into account the finite-volume discretization required by PPM. 
All components of the code are constructed to work with a nonuniform mesh, preserving second- 
order spatial differences. The PPM code uses vacuum boundary conditions for isolated problems, 
preventing inflows when appropriate. The PM code uses a second-order variable-timestep time in- 
tegration scheme. In this paper we have discussed the equations solved by COSMOS and described 
the algorithms used, with emphasis on these features. We have also reported on results from a 
suite of standard test problems, demonstrating that COSMOS works as expected. 

In closing we note some implementation details. COSMOS is designed to run on parallel 
computers using the Parallel Virtual Machine (PVM) message-passing library. A version using 
the Message-Passing Interface (MPI) is currently under development. Also, COSMOS features a 
modular design which simplifies the addition of new physics and the configuration of the code for 
different types of problems. The modular design of COSMOS has served as the basis for FLASH, 
an adaptive mesh-refinement PPM code currently under development at the University of Chicago 
ASCI Flash Center. Details of the FLASH code design will appear in a forthcoming paper (Fryxell 
et al. 1999). 
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Level 1: 
Level 2: 
Level 3: 



6 *~»7-*»* 



x 



Fig. 1. — Example multigrid hierarchy for a 2D grid which is nonuniform in the x-direction and 
uses three levels of grid refinement. The thickness of the lines bounding a cell indicate its grid level; 
for example, the thickest lines indicate cell divisions on the coarsest level. The x-coordinates of cell 
centers on each level are indicated by numbers. 
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Fig. 2. — Schematic diagram of steps in our multigrid solver. Ascending line segments indicate 
prolongation, and descending line segments indicate restriction. Squares indicate direct solution on 
the coarsest level, triangles indicate computation of the residual, and circles indicate correction of 
the residual. 
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Fig. 3. — Density, pressure, velocity, and specific internal energy in the Sod problem at t = 0.206. 
Solid lines represent the analytic solution; squares indicate numerical solution values interpolated 
along the line normal to the shock which passes through the origin. 
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Fig. 4. — Density, pressure, and velocity for the Sedov explosion problem on a uniform 64 3 grid at 
t = 0.0508. Squares represent the angle-averaged numerical solution, and solid lines represent the 
analytic solution. 
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Fig. 5. — Percentage deviation from the uniform-grid case of the angle-averaged solution to the 
Sedov explosion problem on grids of varying nonuniformity. 
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Fig. 6. — Test of the Jeans dispersion relation for gravitationally stable perturbations. Squares 
indicate oscillation frequencies taken from the test calculations; the solid line indicates the expected 
relation. 
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Fig. 7. — One-dimensional convergence results for the gas-only ZePdovich pancake problem, as 
functions of the number of zones N, for fixed perturbation wavelength A. Solid lines indicate the 
LI error norm (eq. [50]), while shaded areas indicate RMS deviations of the error from the LI norm. 
Dashed lines indicate power-law relations for comparison. 




Fig. 8. — Comparison of the 256 2 gas-only pancake test at three different redshifts with the 256- 
zone ID solution. Density, temperature, and velocity amplitude are shown. Symbols represent 
values interpolated from the 2D calculation; lines represent values taken from the ID calculation. 
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Fig. 9. — Average error in Miyamoto-Nagai potential calculations for different flatness ratios b/a, 
plotted cis ct function of the largest multipole moment ^max 

used in the computation of isolated 

boundaries. 
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Fig. 10. — Results of the PM force resolution test, (a) Mean binned radial (filled circles) and 
azimuthal (open circles) components of the interparticle force for randomly positioned particle 
pairs. The dashed line indicates the expected 1/r 2 force. Here r is expressed in units of the zone 
spacing Ar. (b) RMS deviations for the binned radial and azimuthal components, expressed as 
fractions of the mean binned radial force. 
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Fig. 11. — Results of the PM force resolution test with a single pair of particles, where the first 
member of each pair is placed at the center of the grid and the second member is placed at increasing 
distances from the center along the line connecting the center to one of the grid boundaries. The 
plot shows the ratio of the force on the second particle to the expected value as a function of the 
ratio of the particle separation r to the box size L. 
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Fig. 12. — RMS error in the circular orbit of two particles as a function of the ratio of the zone 
spacing Ax to the (fixed) interparticle separation r. Choices of timestep At corresponding to 0.015 
orbital period (triangles) and 0.046 orbital period (squares) are shown. 
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Fig. 13. — Time evolution of the radius of different mass shells in the dark-matter-only spherical 
collapse test. Units are such that the zone spacing is unity. Points are results from a 32 3 simulation; 
solid lines are the analytic solution, plotted until the time of maximum physical radius. 
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Fig. 14. — RMS errors in position (Axrms) and velocity (A-urms) for the dark-matter-only plane 
wave problem in the linear regime. Lines represent results from simulations using a 16 3 grid and 
16 3 particles. Results for three different values of the wavenumber k are shown. Points represent 
values from equivalent runs by Efstathiou et al. (1985), who used an FFT Poisson solver instead of 
the multigrid method used here. 



